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Abstract 

We derive formulas for the Coulomb matrix within the full-potential linearized augmented-plane-wave (FLAPW) 
method. The Coulomb matrix is a central ingredient in implementations of many-body perturbation theory, such as 
the Hartree-Fock and GW approximations for the electronic self-energy or the random-phase approximation for the 
dielectric function. It is represented in the mixed product basis, which combines numerical muffin-tin functions and 
interstitial plane waves constructed from products of FLAPW basis functions. The interstitial plane waves are here 
expanded with the Rayleigh formula. The resulting algorithm is very efficient in terms of both computational cost 
and accuracy and is superior to an implementation with the Fourier transform of the step function. In order to allow 
an analytic treatment of the divergence at k = in reciprocal space, we expand the Coulomb matrix analytically 
around this point without resorting to a projection onto plane waves. Without additional approximations, we then 
apply a basis transformation that diagonalizes the Coulomb matrix and confines the divergence to a single eigenvalue. 
At the same time, response matrices like the dielectric function separate into head, wings, and body with the same 
mathematical properties as in a plane-wave basis. As an illustration we apply the formulas to electron-energy-loss 
spectra (EELS) for nickel at different k vectors including k = 0. The convergence of the spectra towards the result 
at k = is clearly seen. Our all-electron treatment also allows to include transitions from 3s and 3p core states in 
the EELS spectrum that give rise to a shallow peak at high energies and lead to good agreement with experiment. 
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1. Introduction 

For the ab initio calculation of electronic excitations and spectroscopic functions, where variational ground- 
state schemes like Kohn-Sham density-functional theory [l] are not strictly applicable, many-body pertur- 
bation theory has now become the method of choice in applications to solids and their surfaces. It is based 
on a Green-function formalism and an adiabatic switching-on of the Coulomb interaction [2]. In this way 
the Green function of the fully interacting many-electron system can be expanded in powers of the Coulomb 
potential, generating a series of Feynman diagrams with increasing complexity. Practical approximations 
can be designed by terminating the series at a given order or restricting the summation to certain classes 
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of self-energy diagrams that describe dominant scattering processes. A prominent example is the exchange- 
only Hartree-Fock approximation, which includes all electronic interaction effects up to linear order in the 
Coulomb potential, while additional correlation effects resulting from dynamic screening in an itinerant elec- 
tron system are taken into account in the GW approximation [3j. The latter has been successfully applied 
to a variety of materials, especially semiconductors, and generally yields electronic band structures and 
quasiparticle properties in very good agreement with experimental data [4]. The dielectric function, which 
already appears as an intermediate quantity in the GW approximation, can be expanded in a similar manner 
and is itself the key quantity for the theoretical description of optical absorption and related spectroscopies. 

For a numerical evaluation of self-energies or dielectric functions in diagrammatic terms it is necessary to 
project the Coulomb potential, as well as Green functions and all other relevant propagators, onto a suitable 
basis set. Within the diagrammatic expansion the Coulomb interaction describes the elastic scattering of 
two electrons or holes, with a possible momentum transfer between initial and final states. The basis for 
the matrix representation of the Coulomb potential must hence be able to properly describe products of 
initial-state and final-state wave functions. So far most practical implementations have employed a plane- 
wave basis set in combination with norm-conserving pseudopotentials. As the product of two plane waves is 
again a plane wave, this approach has the advantage that products of wave functions can easily be expressed 
in the same basis as the original wave functions themselves. In addition, fast Fourier transformations may 
be exploited, and the Coulomb matrix in reciprocal space is known analytically. For semiconductors, in 
particular, sophisticated theoretical calculations of optical absorption [5j and electron-energy-loss spectra 
[6], which also include excitonic contributions, have been performed in this way. 

While the plane-wave pseudopotential approach works well for sp-bonded semiconductors and simple 
metals, it becomes inefficient for transition metals and rare earths, where a large number of plane waves 
are needed to accurately describe the localized d or / orbitals. A similar problem occurs in oxides and 
other compounds involving first-row elements due to the hard pseudopotentials that only contain minimal 
screening of the ionic core by the innermost Is electrons. Therefore, these materials are best studied within 
an all-electron scheme that treats core and valence shells on an equal footing and already incorporates 
the rapid oscillations of the wave functions close to the nuclei in the basis functions themselves. Here we 
focus on the full-potential linearized augmented-plane-wave (FLAPW) method [7], which is widely used for 
electronic-structure calculations of such materials. It divides space into nonoverlapping muffin-tin spheres 
centered at the atomic positions and into the interstitial region. Inside the muffin-tin spheres the basis 
functions are constructed from numerical solutions of the radial Schrodinger equation with fixed energy 
parameters, whose products lie outside the vector space spanned by the original basis functions. Therefore, 
products of the original basis functions may instead be used to construct a mixed product basis [8], in which 
the matrix elements of the Coulomb potential with initial and final states are then accurately represented. 

While the Coulomb matrix is diagonal in a plane-wave basis and given by a simple analytical expression, its 
evaluation in the mixed product basis of the FLAPW scheme is much more cumbersome. First, the matrix is 
no longer diagonal, and all elements must be calculated numerically. This requires an efficient computational 
procedure. Second, due to the long-range nature of the Coulomb potential v(r) — 1 jr in real space, the matrix 
diverges in the limit of small wave vectors k — > 0. Whereas this divergence is confined to the single head 
element in the case of a plane-wave basis, all matrix elements now contain divergent terms proportional to 
l/k 2 and 1/k. Previous all-electron implementations |9I10| of many-body perturbation theory have often 
bypassed this problem by reverting to a plane- wave basis for the Coulomb potential and related propagators, 
such as the dielectric function, but the projection leads to a loss of accuracy, because the rapid oscillations 
of the orbitals close to the atomic nuclei cannot be resolved then. As a consequence, physical effects like core 
polarization are inadequately described. An alternative approach, the so-called offset-r method, employs an 
auxiliary k-point mesh that is shifted from the origin by a small but finite amount |8|11| . In this way it avoids 
the singularity, but the use of additional meshes increases the numerical cost; even in the most favorable 
case, for cubic symmetry, the number of k points must at least be doubled. Furthermore, the convergence 
of Brillouin-zone (BZ) integrals involving the Coulomb matrix, for example for the GW self-energy, may be 
slow with respect to k-point sampling due to the approximate treatment of the quantitatively important 
region near the zone center. 

In this work we derive formulas for the Coulomb matrix in the mixed product basis including its math- 
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ematically exact expansion around k = 0, which involves terms proportional to 1/fc 2 and 1/fc as well as 
constant terms. A proper treatment of the small- wave- vector limit is especially important for the theoret- 
ical description of optical spectroscopies with zero momentum transfer, but also for the calculation of the 
nonlocal Hartree-Fock potential or the GW self-energy, which both involve an integration over the BZ. In 
a second step, to simplify the numerical treatment we then apply a basis transformation that diagonalizes 
the Coulomb matrix. This eliminates all 1/fc terms and again restricts the 1/fc 2 divergence to a single di- 
agonal element, belonging to a constant eigenfunction. The final situation is thus once more analogous to 
a plane-wave representation, where the dielectric function naturally decomposes into head, wing, and body 
elements, but we retain the full accuracy of the FLAPW basis set. Furthermore, the present algorithm is 
very efficient; the computational time for a well-converged Coulomb matrix with 10 5 elements takes less 
than a second on a modern single-CPU personal computer. The present algorithm is implemented in Spex 
|12j . a computer code for the calculation of excitation spectra and quasiparticle energies within the GW 
approximation. 

This paper is organized as follows. Section [2] shortly describes the FLAPW method and the mixed product 
basis used in this work. The formulas for the Coulomb matrix at finite wave vectors are derived in Section [31 
We then discuss its expansion around k = and the subsequent diagonalization in Section 01 As a practical 
illustration, in Section [5] we present electron-energy-loss spectra of Ni calculated at finite k vectors as well as 
k = within the random-phase approximation. Finally, Section [6] summarizes our main conclusions. Unless 
stated otherwise we use Hartree atomic units. 



2. Basis sets 



2.1. FLAPW method 



In the FLAPW method space is divided into nonoverlapping atom-centered muffin-tin (MT) spheres 
and the interstitial region (IR). The core-electron wave functions, which are (mostly) confined to the MT 
spheres, are directly obtained from a solution of the fully relativistic Dirac equation. The valence-electron 
wave functions with spin a are expanded in interstitial plane waves (IPW) in the interstitial region and 
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for p = and their first energy derivatives u^ lml (r) — du° lm0 (r) / de°i for p = 1, where V eff (r) is the 
spherical average of the effective potential, are suitably chosen energy parameters, and Yi rn (e r ) denote 
the spherical harmonics. The notation e r = r/r with r = |r| indicates the unit vector in the direction of r. 
In a given unit cell the Kohn-Sham wave function at a wave vector k with band index n and spin a is then 
given by 
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with the crystal volume V, the number of unit cells N, and cutoff values Z max and G max . The coefficients 
A-ahnp are determined by the requirement that the wave function is continuous in value and first radial 
derivative at the MT sphere boundaries. If desired, additional local orbitals [14] or higher energy derivatives 
[15] can also be incorporated by allowing p > 2. 
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2.2. Mixed product basis 



The FLAPW method uses continuous basis functions that are defined everywhere in space but have a 
different mathematical representation in the MT spheres and the IR. For the expansion of wave-function 
products, however, it is better to employ two separate sets of functions that are defined only in one of the 
spatial regions and zero in the other. In this way, linear dependences that occur only in one region can 
easily be eliminated, which overall leads to a smaller and more efficient basis. The resulting combined set of 
functions is called the mixed product basis. 

Inside the MT spheres the mixed product basis must accurately describe the products 

l+V L 

Klmp{ r )Kl'm'p'( r ) =Kl P ( r ) Y hn( e r)<l>p'( r ) Y l'rn>(e r ) = ^ ^ C lml'm> LmK LP (r)Y LM (e r ) , (3) 
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which we expand in spherical harmonics with the Gaunt coefficients 

Cimi'm'LM = J Y{^ (e r )Y^ ro / {g t )YI m (e r ) dQ . (4) 

The index P counts the radial product functions U° LP (r) = Uaip( r ) u ai' p' (r) for a given angular quantum 
number L. We emphasize again that, in general, the latter lie outside the vector space spanned by the original 
numerical basis functions {u^ lmp (r)}. Initially, the set of radial product functions is neither normalized nor 
orthogonal and usually has a high degree of (near) linear dependence. An effective procedure to remove these 
(near) linear dependences is to diagonalize the overlap matrix and to retain only those eigenvectors whose 
eigenvalues exceed a specified threshold value [13J. In this way the MT functions become orthonormalized. By 
using both spin-up and spin-down products in the construction of the overlap matrix we make the resulting 
basis spin-independent. If desired, the basis set may be reduced further by introducing an additional cutoff 
value L max for the angular quantum number. On the other hand, it must be supplemented with a constant 
MT function for each atom in the unit cell, which is later needed to represent the eigenfunction that 
corresponds to the divergent eigenvalue of the Coulomb matrix in the limit k — > 0. From the resulting 
orthonormal MT functions M a LUp{ T ) — M a Lp(r)YLM i^r) we formally construct Bloch functions 

W = ^E^'^^^^^-T-^D^^-T-Ro)- (5) 

The sum runs over all lattice translation vectors T, and M a i,p(r) = if r is larger than the muffin-tin radius 

Sa- 
in the IR, since the product of two IPWs equals another IPW, we use the set 

M G (r) = J= e *( k + G )-e(r) (6) 
if r € MT 

6(r) = < (7) 
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and a cutoff G' mgx < 2G max in reciprocal space. Together with the MT functions we thus obtain the mixed 
product basis {Mj (r)} = {M a \ M p(r), Mg(r)} for the representation of wave-function products. In contrast 
to the MT functions, which were explicitly orthonormalized, the IPWs are not orthogonal to each other; the 
elements of their overlap matrix can be calculated analytically and are given by 

(m£|m&) = ^o GG -(k) = <we G _ G , , (8) 

where 
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are the Fourier coefficients of the step function ([7]) and fl denotes the unit-cell volume. We also define a 
second set, the biorthogonal set 

M / k (r)=^07 / 1 (k)M,^(r) (10) 
J 

with the overlap matrix Ojj(k) = (Mf \Mj). It fulfills the identities 
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where the completeness relation is only valid in the subspace spanned by the mixed product basis, however. 
As the MT functions and the IPWs are defined in different regions of space and the MT functions are 
orthonormal, only the IPWs overlap in a nontrivial way. It should be noted that the overlap matrix is 
k-dependent because the size of the mixed product basis varies for different k vectors. 

For the evaluation of the Coulomb matrix elements we have to find a numerically tractable expression for 
the IPWs. A straightforward approach might employ the Fourier transform of the step function and rewrite 
H} as a sum over reciprocal lattice vectors 

h(r)= lim ^= E e G ,e*( k+G+G '> r , (12) 
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where Gpw is a cutoff radius in reciprocal space, for which we must of course choose a finite value in practice. 
Eq. l[T2"Tl has a very simple mathematical structure and is easy to implement. For example, the calculation of 
the matrix elements IPW-IPW only involves Fourier coefficients of the step function 6(r) and the Coulomb 
interaction 1/r, which are both known analytically. As an alternative, we may exploit the Rayleigh expansion 
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involving the spherical Bessel functions ji (x) in order to subtract the plane waves inside the MT spheres 
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where we use the abbreviations q = k + G and r' = r T R a , and 9(r) denotes the Heaviside function. 
In a practical implementation we must use a finite maximal angular momentum (pw, which thus becomes 
the relevant convergence parameter. Despite its more complicated mathematical appearance, we have found 
that this representation in fact facilitates a considerably faster numerical evaluation because of the slow 
convergence of the step function in f(T2|) with respect to the number of Fourier coefficients. We illustrate this 
point in Section [5l In our subsequent derivation we hence employ expression (fl4| . 



3. Coulomb matrix at finite k 



In this section we derive the formulas for the computation of the Coulomb matrix elements 

viAk)= jjMrmm dW (15) 

for arbitrary finite wave vectors; the limit k — > is discussed in Section 21 Due to the composite basis 
set {Mf(r)}, which consists of MT functions with / = (aLM P) and IPWs with / = G, the Coulomb 
matrix is made of four distinct blocks. As it is Hermitian, however, the two off-diagonal blocks are complex 
conjugates of each other, and thus we have to consider only three blocks explicitly, which correspond to the 
combinations MT-MT, MT-IPW, and IPW-IPW. Svane and Andersen [16|| already examined the matrix 
elements MT-MT for finite k vectors in the context of the linearized muffin-tin orbital (LMTO) method. In 
the following we summarize the derivation in a somewhat different notation and then give the expressions 
for the additional matrix elements involving IPWs. 
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3.1. MT-MT 



If we insert the Bloch representation ((5j) for the MT functions in 
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where the difference vector R aa ' = R a ' — R a points from one MT center to another in the same unit cell. 
The integrals in (fl7|) corresponding to R = and R^O with R = T + R aQ < give rise to two contributions 

VaLMP,a>L'M>P>(k) = S aa 'V^LM P,aL> M> P' + ^aLMP^'L'M'P'O^ ' ( 18 ) 

which we evaluate separately in the following. 

Let us first consider the integral for R = 0. It can be simplified considerably by inserting the identity 
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where r< and r> indicate the smaller and larger value of {r, r'}, respectively. After carrying out the angular 
integrations we obtain 
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The remaining integrations can be easily performed by standard numerical techniques on a radial mesh. 
For the integrals with R^Owe may formally define a multipole potential 
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that acts in the first MT sphere as a result of a "charge distribution" M a > l'm> p> (r' — R) in the second, where 
Qa> l> p> denotes the multipole moments 
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the multipole potential (|2lj) created by a MT function at R can then be written in terms of radial functions 
and spherical harmonics at the origin. The corresponding "electrostatic interaction energy" is given by 
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After performing the sum over lattice vectors in ([17]) we eventually obtain 
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with 

St (k) = ]T ^ T ,„ T * |t+i YL(e T+Raal ) , (27) 

where the sum runs over all lattice vectors excluding the case T + R aa / = 0. We note that Sf^ (k) is 
closely related to the structure constant defined in the context of the LMTO method [17]; however, it is 
not dimensionless and therefore not a constant of a given crystal structure. For the numerical evaluation of 
SjJjjJ (k) one must apply the Ewald summation technique. 

3.2. MT-IPW 

For the matrix elements in the off-diagonal block 

WBOM - // ^ff" *r*S (28) 
we can again introduce a formal "charge distribution" given by Mg(r') that creates a potential 

$( r )= / f^idV = 4= 4 ^ 2-/ n rr^V , (29) 
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where the integral runs over the combined volume of all MT spheres, cutting out the plane waves inside 
them. The "electrostatic interaction energy" arising from the first term in the brackets is given by 

■44 / M Q V M p(rK qr d 3 r = il^F L V(e q )e jGR »l f" r*M aLP {r)j L (qr) dr , (30) 
V V 1 J v il 9 Jo 

where we have again used the Rayleigh expansion (fT3"ll and the abbreviation q = k + G. If the exponential 
function in the second term on the right-hand side of f29|) is also replaced by the Rayleigh expansion, then 
the resulting integrals are equivalent to those considered in Section I3TT1 above. We can hence evaluate them 
in the same way. The resulting final expression for the Coulomb matrix element 

VaLMP, G (k) = 4iMF,G( k ) + V aLMP,G ( k ) + V aLMP,G ( k ) ( 31 ) 

consists of three distinct terms, which are given by 
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with the multipole moments 
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and the integrals 

Hl,r) = j r r> l+2 3 l (.qr>)dT> and J al (q,r) = J'' dr> , (34) 

for which analytic expressions are given in appendix [Bj 

3.3. IPW-IPW 

The remaining integrals 

voa/(k)= jjHKi^n M , (35) 

are evaluated in a similar manner. The subtraction of the plane waves inside the MT spheres now leads to 
a decomposition of the matrix elements into four terms 

v GG ' (k) = vg G , (k) - „gk (k) - vg G , (k) + wg>,, (k) . (36) 
The first three can be calculated analytically and yield 
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by again replacing the exponential functions with the Rayleigh expansion (fl~3| and following the procedure 
outlined in Section |3~T1 above. The subsequent summation over MT spheres and angular momenta yields 
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with q = k + G, q' = k + G' and the double integral 



Kal(q,q') = £ £ ' r>r*j l (qr)j l (q'r')^ T drdr' . (39) 
For the latter an analytic formula is derived in appendix iBl 

4. Expansion around k = 

Due to the long-range nature of the Coulomb interaction v(r) = 1/r in real space, its Fourier transform 
4ir/k 2 diverges for k — * 0. As a consequence, the Coulomb matrix in the mixed product basis also diverges 
with a leading term proportional to 1/fc 2 . Since the MT functions contain nontrivial k-dependent coefficients, 
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we further have additional terms proportional to 1/k, It is helpful to identify all relevant terms in advance. 
For this purpose we formally represent the basis functions by their Fourier transforms 

M k (r) = 4r c/ G (k)e*( k + G )- r (40) 

with the coefficients 



c /G (k) = -L / e- l ( k+G ) r A/ / k (r)d 3 r. (41) 
v V J 

The sum runs over all reciprocal lattice vectors G. For the IPWs the coefficients are k-independent and equal 
CG'G(k) = 6g-g' for Mg,(r), but for the MT functions they exhibit a nontrivial k dependence. Using the 
expansion 

M k (r) ~ -L £ (c IG + k • Vc IG + ik T Ac /G k) e ^ +G > r (42) 

for k — > with c/g = cr G (k)| k=0 , Vc/ G = VkC/G(k)| k=0 and Ac/ G = VkV k c/G(k)| k=Q , we can write 
the Coulomb matrix elements in this limit as 

47T 47T 

wjj(k) ~ c /o c ^otj + [cj ( e k • Vcjo) + (e k • Vc* I0 )cj ] — + [(e k ■ VcJ ) (e k • Vc 70 ) 



-^Cj (e£Ac 70 e k ) + -(e£Acj e k )cj 



4tt+ 5>J G c JG — -J. (43) 

G#0 M 



Evidently, all matrix elements contain divergent contributions proportional to 1/k 2 in addition to a constant 
term. Furthermore, if Cjc(k) or cjG(k) are truly k-dependent, i.e., for matrix elements that involve MT 
functions, we also have terms proportional to Y 1 * m (ek)/fc and 5^n( ek ) arising from the first and second 
square bracket, respectively (compare (|A.5[) and (| A.8[> ) . As a consequence, we can write 

(44) 

2=0 m =-l 

and from (fl~5|l follows 

«S?=# and «« m = • (45) 

We will see in Section 14.41 below that the terms corresponding to I > can in fact be eliminated if we 
perform a basis transformation to the set of eigenvectors of the Coulomb matrix. Nevertheless, for the sake 
of completeness we will here give the appropriate formulas for Vj)i m with I > in the original mixed 
product basis as well. As in the previous section, we proceed by discussing the blocks MT-MT, MT-IPW 
and IPW-IPW individually. 



4.1. MT-MT 



The second term on the right-hand side of lfT8]). explicitly given in (J26j) , diverges for k — > and L + L' < 2, 
because the leading term of (k) is proportional to k 2 ~' , which is seen in the following way: For small 
k the sum over T in (|27|l is dominated by contributions belonging to large lattice vectors. Then one can 
approximate the sum by an integral 

i f pik-T a -i 

S^(k)~ e -*- R «'- J ^Ctex) d 3 T = - _ e - <k - R "'I£ t (e k )fc f - a , (46) 

where we have used (fT3|) , (|B.2|) , and (| A.4[) . The same expression appears in the first term of the reciprocal- 
space sum corresponding to G = in the Ewald summation. The remaining terms and the real-space sum 
yield an additional constant term Sf^ , so that we obtain 

sz (k) ~ (2i-i)i\n e ~ lk ' R " a ' Yl * Mkl ~ 2 + sz (47) 
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for I <2. After inserting this expansion into l(26|) we obtain 



X „,( a ) I I i\L'+M' n n nan' 

V aLMP,a'L'M'P' —°a,a'V aLM p a , L , M ,p, +\ 1 ) Cl> M> ,LMH a> V P> W aLP o (L+L')(M- W) ' 



u aLMP,a'L'M'P 

with I < 2. 



■ l m — &l,L+L'&m,M-M' (^1) L +M 



Am 1 



(2/-l)!!0 



CL'M'.LAtQa'L'P'QaLP 



(48) 
(49) 



4.2. MT-IPW 



The case MT-IPW is more complicated, because higher orders in the multipole moments QJii m i must be 
taken into account when multiplying with the divergent S?£ +l ,^, M _ m i\ (k) in (|32c|) . In particular, we need 
the expansions of G^Jm (QaAm 1 ) U P to second (first) order 



?*o + o G ~V4^ 



O k + G ~ 47ri— 



ji{Gs a ) - j2{Gs a )s a ke k ■ e G + ^(Gs^s 2 ,*; 2 (e k ■ e G 



,2 lj 2 (Gs a ) s fc2 



r im( e G) (j2(Gs Q ) - j 3 (Gs a )s a fce k • e G ) + Y* m (e k ) 



G 



2 G 



For G = these simplify to 
\f\~n 



Qa 



) k 

' aim 



00 g 3 a 



3 'i--d 2 



i 

itT a ' 
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(50a) 
(50b) 

(51a) 
(51b) 



Here we have used the identities (|A.1[1 (|A.7[) . In addition, (|32ap contributes to v aLMPG if G = 0. The final 
expression for i)fj MP G and v^ MP G is written as 



(0) _ (0a) (0b) 

V aLMP.G ~ V aLMP,G V aLMP,G ' 



where the quantity 

J 0a ) 



aLMP.G - V 1 - 6 G0)v a i M P,G(°) + ^LMP, G (°) 



,( b ) 



(52) 
(53) 



^| "Jr VV ° 

/^QaLP^2 X! + " 1 ' R "' c l'm',LMQ3l'm' S (L+l' 



)(M-m') 



i'=0 7 



is directly obtained after replacing Sf^ (k) by the terms of zeroth order in the expansion (J4TJ) . The second 
term v^mpg results from multiplying the divergent terms in lj47|) with the higher orders of (|50"f in l|32cp 
as well as, in the case G = 0, the term 1/q 2 with the higher orders of ji(qr) in (|32a,[) . After some algebra we 
obtain 

r ( 4 -) 5/2 ^ ^ 



,(0b) 

aLMP.G 



= < 



30^3/2 V 

(4tt) 5 / 2 



ft 3 / 2 



n JG-tl., Sa ' 



s 3 , 

a' a 


\HGs 


G 


2G 


(4tt) 3 / 2 p 


6Vfi Jo 


s 3 , 

a , s a' 


\h(Gs 


G 


2G 



G 



r 4 M a0P (r)dr 



if G ^ and L = , 

if G = and L = , 

if G f and L = 1 , 
otherwise , 



(54) 
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as well as 



(i) 

aLMP,G;lm 



SliS 



U°Mm 



(4tt) 



5/2-L 



(2L+ l)!!ft 3 / 2 



laLP 



00 



JGO 



(47T) 2 Z L Q aLP 

(2L+ l)!!^ 1 ^ 



(55) 



4.3. IPW-IPW 



Finally, in the calculation of the IPW-IPW matrix elements we can use the fact that the square brackets 
in lj43]) vanish. This simplifies the derivation considerably, because all angular-dependent contributions can 



be discarded from the outset, and hence we have v, 



(i) 

GG'Jr. 



for I > 0. We again write 



„(0) _ w (0a) , (Ob) 
GG' — GG' ^ GG' ' 



where the first contribution v GG , is given by the nondivergent terms in lj36]) after replacing Sf^ (k) by 
which yields 



(0a 



!, = (!- 6go) [4 a G '(°) - 4U )} - (1 - S G 'o)v^ G ,(0) 



l^ e i(G'-G).R a g^|^ y;m(eG)F4(eG()/CQi(g)?0 



(56) 



/=0 m 



E ^ ^ ^ ( — 1)' +m £ e ' GR ° e ' G Ra ' C ''m',im ( 9Sm ( 9°;'m"5'(i+i')(m-m') 



I — m— — I V — m' — — / 



a, a' 



Further, by inserting the expansions ([50)1 and ([51]) as well as the k-dependent terms of |47j) into ([38)1 one 
obtains another constant contribution 

( (Ait) 3 - 2 - 2 



„(0b) 
GG' 



o 2 



^ e -iG-R. e iG'.R a ,^ |-^ 2 (G So )i 2 (G' Sa 0SaSa'(e G -eG0 - |ji (Gs a )j 3 (G'*a')4 

if G ^ and G' ^ 0, 



1 ■ tn \ - tn> \ 2 _, h{Gs a )h{G' s a -)s a , h(Gs a )ji(G' 's a >)s a 
-gj3(Gs Q )ji(G s Q /)s Q + — 



2G" 



n 2 

90ft 2 



^ e4G ' R.f|Ii || Jl(G ' Sa ,) _ l i3 (GV)4 

a, a' 

„3 „2 f 2 i | 

E-iG R / b a b a' ) b a ■ / ri \ -In \ 2 , 1 

6 ^l30 Jl(GSa,) ~T8 J3(GSa ' )v + 6 

a, a' 

£-24 (4 + 4) 



2G 
l j 2 (GV) 
6 G> ' 

Ija (<?*„') 



G 



if G = and G' ^ 0, 
if G ^ and G' = 0, 
if G = G' = . 

(57) 

For the calculation of v GG , no we mus t take the divergent terms of |37j) into account and eventually obtain 



a, a' 



V, 



GG',00 = £ ^ GRaelG ' Ra 'QaolQa'oO + [©G-G^GO + <5 G <o) - S G0 S G , ] . (58) 



4.4. Diagonalization 



In a pure plane- wave representation response matrices and similar quantities decompose into head xoo, 
wings xgo, Xog', and body Xgg' with G,G' ^ 0. These behave differently for k — > 0. For the density 
response function, as an example, head and wing elements are quadratic and linear in k, respectively, while 
the body elements remain finite but still exhibit an angular k dependence. As the mixed product basis is 
related to the set of plane waves by means of a basis transformation, these matrix elements will now mix 
in a complicated manner. It is hence desirable to make another transformation that restores the convenient 
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mathematical properties of the plane- wave basis. For this purpose the basis must include a constant function, 
which corresponds to the limit of e lk r / VV for k — ► and is responsible for the decomposition into head, 
wings, and body. Such a basis is given by the functions 

£ k (r) = £i*M k (r), (59) 

where E^j is the Ith component of the /j,th eigenvector of vjj (k) . In this basis the Coulomb matrix (k) 
becomes diagonal, which is also advantageous in matrix multiplications involving i> M „(k), e.g., for the calcu- 
lation of the dielectric function. Furthermore, the eigenvalues are a direct measure for the probability of the 
elastic scattering between two particles. Small eigenvalues thus identify less important scattering processes 
that might be neglected, leading to a smaller and optimized basis set after removal of the corresponding 
eigenvectors. 

Because of the divergent terms in l(44|) the diagonalization in the limit k — > is not trivial. The first 
eigenvector E k corresponding to the divergent eigenvalue ui(k) = Air/k 2 is, however, easy to obtain from 
the analytic projection of e ik ' r /VV on the biorthogonal mixed-product-basis functions 



cj (k) = — -y L * M (e k ) / r 2 j L (kr)M aLP (r) dr for / = (aLMP) , 
Jo 



5gq for I = G . 

(60) 



which in the limit k — > becomes 



y/4Trsl/(3fl) for I = (aOOl) , 
£ w = { S G0 for I = G , (61) 

otherwise . 

Here we have assumed that the constant MT function of atom a is normalized and identified with the index 
(aOOl). The other eigenvectors E° and eigenvalues v M (0) for [i > 1 are obtained by diagonalizing the last 
term of the formal expansion (|43]l 

vu = c *ig,cjg%, (62) 

G^O 

which is unknown so far. The matrix vfj, calculated in the previous section, contains v%j but also the 



spherical average of the second square bracket in f43|) . If we denote this spherical average by wu, then we 
can write 

vu = vfj - wu . (63) 
In order to evaluate wu we introduce the natural basis 



^-(fca; iky), k\ — ^-y 

which allows us to write the k-dependent terms in the second bracket of (|43|) in terms of spherical harmonics 
according to 



k-i = — (k x - iky), k\ = —7={— k x — ik y ), k — k z , (64a) 

V2 

d-i = (d x + id y ) , di = — = {-d x + id y ) , d = d z , (64b) 
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e£Ac J0 e k = E y im( e k)yim'(ek)5^a m -cjo 



(65) 



m— — 1 m'= — 1 
1 



q X] ^m^m c JO + ^ ^lmam2(m-m')*2*(m-m') ( e k)<9„<9 m ' Cj , 



4tt 



(e k • Vc^ ) (e k • Vcjo) = — ^ ^ r i * m (e k )y lm /(e k )S^cJ 9 m /c J0 



(66) 



m— — 1 m'=— 1 

1 



C,70 , 



m=— 1 m=— lm'= — 1 

where we define d m ci = <9 m cj (k)| k=0 and similar abbreviations. The last equation follows from the identity 



e k ■ Vc/o = \/Y £ yim( e k)9 m c/ ■ 

m— — 1 

When we take the spherical average, the harmonics with Z > vanish, and we obtain 

1 



4tt 



wij = 



m— — 1 



(#mCj ) (a m c J0 ) + ^Cj a^9 m Cj + ^C J0 9 m 9^Cj 



(67) 



(68) 



with 



C/O = < 



d m cio - 



CfO 



V47rs|7(3fi) for 7 = (aOOl) , 

e_ G for 7 = G , 

otherwise , 
/~4 f Sa 

-6m5 M m]j J r 3 M alP (r) dr for 7 - (alMP) 

otherwise , 

-5 L o\l / r 4 A7 a0P (r) dr for 7 = (aOMP) , 



(69) 



(70) 



(71) 



otherwise . 



5. Test calculations 



Apart from the evaluation of (|27| . which is easily converged to high precision by means of the Ewald 
summation technique, and the radial meshes for numerical integration inside the MT spheres, the cutoff 
value Ipw is the only convergence parameter in the construction of the Coulomb matrix elements presented 
above. On the other hand, in an alternative implementation that uses the representation lfl2|) rather than 
the Rayleigh expansion for the IPWs the matrix elements must be converged with respect to the reciprocal 
cutoff radius Gpw- FigureQ] compares the convergence behavior of these two approaches. The curves indicate 
the root mean square deviation of the Coulomb matrix for bulk silicon, averaged over all matrix elements 
MT-IPW and IPW-IPW and over 64 k points, from the fully converged results calculated with Zpw = 26. In 
both cases we employ the same cutoff parameters G' m ^ = 3.6Bohr _1 and 7 max = 4 for the mixed product 
basis, the MT functions are constructed from products u^ lQ (r)u^ llQ (r) with I < 2 and I' < 3. On average 
this yields 411 basis functions per k point. It is evident that the results obtained with the present method 
converge much faster than those obtained with the Fourier transform of the step function. Furthermore, at 
the same level of accuracy the present method is typically by a factor of 10-100 faster. 
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As an application we now consider the simulation of experimental spectroscopies related to the complex 
dielectric function, which describes many-body screening effects in a correlated electron system. In electron- 
energy- loss spectroscopy (EELS), for example, the measured differential scattering cross section is directly 
proportional to the imaginary part of a diagonal element of the inverse dielectric function [19] 



V 



s-\r,r';u)e ik < T '- I U 3 rd 3 r' , 



whereas in optical absorption one measures the imaginary part of [20] 

£ M M = lira l/e^(k,u) . 

k— >0 

In the framework of many-body perturbation theory the dielectric function is written as 

e(r, r'; w) = S(r - r') - J v(r, r")P(r", r'; u) d 3 r" 

with the polarization function P(r, r';u>) and the Coulomb interaction v(r,r') = l/|r — r'|. We use the 
random-phase approximation 



(72) 
(73) 
(74) 



occ unocc 



P(r,r';o,) = ££ E ^k(r)^ q+k (r)^k(r')^ q+k (r') 

o n,q n'.k 



(75) 



1 



UJ + e". - e 



n'q+k 



+ ir ) w - e° + e 



n'q+k 



IT] 



where 77 is a positive infinitesimal. As P(r, r';ui) contains products of wave functions evaluated at r and r', 
it can be represented in the mixed product basis as 



P(r,r';u,) = V f P/j(k, uj)Mf(r)Mj* (r') d 3 k 
i,j Jbz 



with complex coefficients 



P(r, r';u))Mf (r)M^(r') d 3 r d 3 r' 



(76) 



(77) 



Next we transform this matrix to the basis given by ([59]) . This yields P M „ (k, lj) , which in the limit k — > 
decomposes into head, wing, and body elements as discussed in Section l4~4l We use the tetrahedron method 
for integrations over the BZ. 



1 tr- 



ier 



1 



2 10"' 

5 10-" 



10" 



10" 1 



10" 



(a) 



x 

< 



10 12 14 16 18 20 

I PW 




10 12 14 16 18 20 
G PW [Bohr" 1 ] 



Figure 1. Average root mean square deviation Av from the converged matrix elements (MT-IPW and IPW-IPW) as functions 
of (a) the convergence parameters ipw and (b) the reciprocal cutoff radius Gpw for the Fourier transform of the step function 
in (1121) . The mixed product basis was optimized for Si bulk. 
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^ 1 1 1 1 1 

20 40 60 80 100 

co[eV] 

Figure 2. EELS spectra of spin-polarized Ni for k = 27r/aNi(£,£,0 with £ = (solid line), £ = 0.1 (dashed line), and £ = 0.2 
(dotted line). 

In the long- wave- length limit k — > we must carefully expand the polarization function around k = 0, 
since it is multiplied with v(r,r') in J74|), which diverges in this limit. Because of the orthogonality of the 
wave functions the projection (£i <^q|<^' q+ k) = ( e * k '"Vnql^n'q+k)/^^ i s linear in lowest order in k for 
interband transitions with n ^ n' . We calculate this leading term with k • p perturbation theory [21]. For 
a metallic system the sum in f75|) also contains contributions from intraband transitions with n = n! at 
k = 0. It can be shown that these are nonzero only for the head element and given analytically by the Drude 
formula [22], which is quadratic in k. The latter depends on the plasma frequency, which we obtain by an 
integration over the Fermi surface. In conclusion, the head and wing elements of P^ u (k,uj) are quadratic 
and linear in k, respectively. If we use the symmetrized dielectric matrix 

i^ik, uj) = «v - «y 2 (k)p^(k, ^) v y 2 (k) , (78) 

1 /2 

where the w M (k) are the eigenvalues of vu(k), then all elements of ^(k, iS) are finite because v x ' (k) = 
\/47r/fc. We note that the diagonal quantities considered above remain unchanged with this symmetrized 
definition. As the first eigenvector of u/j(k) corresponds to the projection of e lk ' T /y/V onto the biorthogonal 
mixed product basis, the head element (k, u) directly equals the spectroscopic function (72]). 

In figure [2] we show EELS spectra Ime _1 (k,w) for spin-polarized Ni at three k vectors 27r/aNi(£, £, £) 
with £ = 0.0, 0.1, 0.2 and the lattice constant aNi = 6.66 Bohr We use the parameters i max = 8, G max — 
3.57 Bohr" 1 for the FLAPW and L max = 4, G max = 5.00 Bohr" 1 for the mixed product basis. The BZ is 
sampled by 1661 points in its irreducible wedge, corresponding to a 40x40x40 k-point mesh in the full zone. 
As the spectrum extends over a wide energy range up to 100 eV, we augment the FLAPW basis by second 
and third energy derivatives as local orbitals to guarantee an accurate description of high-lying conduction 
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T 




co[eV] 

Figure 3. EELS spectra of spin-polarized Ni for k = 27r/aNi(0.25, 0, 0) with (solid line) and without (dashed line) core-state 
contributions. The inclusion of transitions from core into conduction states gives rise to a shallow peak starting at 62 eV, which 
is also seen in experiment (symbols) |23j . 

states [H]. For the construction of the MT functions M a Lp(r) we employ products uZi P ( r ) u ai' P ( r ) with 
I < 2, I' < 3, and p — 0, i.e., energy derivatives (p > 1) are neglected. In the calculation of ([75)1 we take 118 
conduction and the 10 valence states as well as the eight 3s and 3p core states into account. As we invert 
the dielectric function, local-field effects are fully included. As seen from the figure, the spectra are very 
similar for the three k vectors. When compared with the curves calculated at finite k points, the spectrum 
for k = clearly constitutes the limit k — > 0. 

As already pointed out, the spectra in figure [2] already include transitions from the 3s and 3p core states 
into conduction states. Figure [3] shows a comparison of spectra calculated with (solid line) and without 
(dashed line) these core-state contributions at k = 27r/aNi(0.25, 0, 0). The largest difference between the 
two curves is seen around 62 eV, which corresponds roughly to the threshold energy required to excite a 
3p electron above the Fermi level. These additional transitions give rise to a shallow peak, which is also 
observed in experiments with an onset at the same energy. The inclusion of transitions from 3s and 3p core 
states into conduction states within our all-electron method thus brings the calculated spectrum very close 
to experiment (symbols) |23j . 

6. Summary 

In this work we have derived formulas for the Coulomb matrix elements within the all-electron FLAPW 
method. As the Coulomb interaction couples two incoming and two outgoing states, a suitable basis set 
must be capable of accurately represent wave-function products. Such a set is given by the mixed product 
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basis, which contains MT functions as well as interstitial plane waves. We use the Rayleigh expansion for the 
latter, because it makes a very efficient numerical implementation possible. Furthermore, we have derived 
an exact expansion of the Coulomb matrix around k = that isolates all divergent terms ~ k~ 2 and ~ fc _1 . 
Most of these vanish if we then make a basis transformation to the eigenvectors of the Coulomb matrix. The 
properties of this new basis set are formally similar to those of a plane-wave basis. In particular, response 
functions decompose into head, wing, and body elements with the same characteristic dependence on k. 
However, the basis construction of this involves no approximation, and the accuracy of the FLAPW basis 
set is completely preserved. 

As an illustration we have shown EELS spectra for ferromagnetic Ni at different k vectors including k = 
0. Very good agreement with experiment was achieved over a large energy window by taking core-electron 
contributions into account. 
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Appendix A. Mathematical relations 

In the derivations of Sections [3] and [4] we have used the following relations: 

ji-i(x) + j t+1 {x) = (21 + l)ji(x)/x , 



d 

dx 
d 

dx 



3i(x) = -ji{x) - j l+ i(x) , 

X 

l + l 



M x ) =3l-i[x) 

J. 

Ji{x) = 



(2Z + 1)!! 



1 - 



4Z + 6 



+ 0(x 4 ) , 



e k -e G = — Y y im( e k)i / i* m (eG), 



*l*m( e k+G) = Y* m (e G ) 



Cljiiljn'2(m'-m)^2,m'-m(^a) 3 
1 



2k 



3G 



m f — — l 



(e k • e G ) 2 = i + ^ J2 Y 2 * m ,(e k )Y 2m ,(e G ), 



Y lm (e k ) = \ — 



3 15 

3 fart 

47T k 



m' = -2 



1 1 

J2 d m dU(k)Y lm (e k ) = pyim(ek) [dk (k 2 d k ) - 1(1 + 1)] /(A;) . 



o(fc 2 ) 



A.l 
A. 2 

A.3 

A.4 

A.5 

A.6 

A.7 

A.8 
A.9 

(A.10 



Appendix B. Integrals over spherical Bessel functions 

The derivations in Section [3] give rise to a number of integrals over spherical Bessel functions that can be 
evaluated analytically. Explicit formulas for ([3~4"ll follow from the recursion relations l|A.l|) - (|A.3j) 
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J+2 



-ji+i(qr) if q^O, 



&nrr 3 



if 9 = 0, 



Jai{q,r) 




j=T3l-i(qr) - -j^ji-iiqsa) 



Sio^(s 2 a -r 2 ) 



if q = . 



(B.l) 



(B.2) 



We can also find an analytic expression for the double integral f39|) , because the above integration formulas 
and the recursion relation l|A.ip lead to the solution 



r t '\ 21 + 1 
Kal{q,q) = ^- 



q 2 _ q ,2 



r 2 ji{qr)ji(q'r) dr ^ji+i(qs a )jl-i(q's a ) 

qq 

—ji+i{qsa)ji-i{q'sa) — jji-i(qsa)ji+i(q's a ) 
q q 



3i+i{qsa)ji+i(q's a ) 21 + 1 ji +2 (qs a )ji(q's a ) - ji(qs a )ji +2 (q' s a ) 



qq' 



21 + 3 



2 O 



(B.3a) 
(B.3b) 



where we used the symmetry of K, a i(q,q') with respect to q and q' to eliminate Jj° r 2 ji(qr)ji(q'r)dr. The 
expressions l|B.3aP and (|B.3b|l are stable for large and small q, q' , respectively. The limiting cases are 



lim IC a i(q,q') = Si a -^ [qs a ji(qs a ) + h{qs a )] for q ^ 



lim )C al (q,q') = [(21 + 3)jf +1 (qs a ) - (21 + l)ji(qs a )j l+2 (qs a )] for q± , 

q'^q 2q l 

2 

lim Kai (q, q) = 5 m — sl . 

q,q'^0 15 



(B.4a) 
(B.4b) 
(B.4c) 
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